﻿using RayDen.Library.Core;
using RayDen.Library.Core.Primitives;

namespace RayDen.Library.Components.Volume {
    public static class PhaseFunctions {
        public static float RdIntegral(float alphap, float A) {
            float sqrtTerm = MathLab.Sqrt(3f * (1f - alphap));
            return alphap / 2f * (1f + MathLab.Exp(-4f / 3f * A * sqrtTerm)) *
                   MathLab.Exp(-sqrtTerm);
        }


        public static float RdToAlphap(float reflectance, float A) {
            float alphaLow = 0f, alphaHigh = 1f;
            float kd0 = RdIntegral(alphaLow, A);
            float kd1 = RdIntegral(alphaHigh, A);
            for (int i = 0; i < 16; ++i) {
                float alphaMid = (alphaLow + alphaHigh) * 0.5f;
                float kd = RdIntegral(alphaMid, A);
                if (kd < reflectance) { alphaLow = alphaMid; kd0 = kd; }
                else { alphaHigh = alphaMid; kd1 = kd; }
            }
            return (alphaLow + alphaHigh) * 0.5f;
        }
        public static float PhaseIsotropic(ref Vector w, ref Vector wp) {
            return 1f / (4f * MathLab.M_PI);
        }

        public static float PhaseRayleigh(ref Vector w, ref Vector wp) {
            float costheta = Vector.Dot(ref w, ref  wp);
            return 3f / (16f * MathLab.M_PI) * (1f + costheta * costheta);
        }

        public static float PhaseMieHazy(ref Vector w, ref Vector wp) {
            float costheta = Vector.Dot(ref w, ref wp);
            return (float)((0.5f + 4.5f * MathLab.Pow(0.5 * (1f + costheta), 8f)) / (4.0 * MathLab.M_PI));
        }
        public static float PhaseMieMurky(ref Vector w, ref Vector wp) {
            float costheta = Vector.Dot(ref w, ref wp);
            return (float)((0.5f + 16.5f * MathLab.Pow(0.5 * (1f + costheta), 32f)) / (4.0 * MathLab.M_PI));
        }
        public static float PhaseHG(ref Vector w, ref Vector wp, float g) {
            float costheta = Vector.Dot(ref w, ref wp);
            return 1f / (4f * MathLab.M_PI) *
                   (1f - g * g) / (float)MathLab.Pow(1.0 + g * g - 2d * g * costheta, 1.5f);
        }
        public static float PhaseSchlick(ref Vector w, ref Vector wp, float g) {
            float alpha = 1.5f;
            float k = alpha * g - (1f - alpha) * g * g * g;
            float kcostheta = k * Vector.Dot(ref w, ref wp);
            return 1f / (4f * MathLab.M_PI) *
                   (1f - k * k) / ((1f - kcostheta) * (1f - kcostheta));
        }
    }
}